/*==============================================================================
案例1：员工薪酬影响因素分析
文件: case01_salary_analysis.do
==============================================================================*/

clear all
set more off
capture log close

* 设置工作目录
cd "/Users/mac/git/stata"

* 创建输出目录
capture mkdir output
capture mkdir output/figures

* 开始日志
log using "output/case01_salary_analysis.log", replace text

display "=========================================="
display "案例1：员工薪酬影响因素分析"
display "=========================================="

/*------------------------------------------------------------------------------
第一步：数据加载和探索
------------------------------------------------------------------------------*/

display ""
display "第一步：数据加载和探索"
display "----------------------"
display ""
display "📊 知识点："
display "  - 数据加载：从网络或本地加载数据"
display "  - 数据清理：处理缺失值"
display "  - 变量转换：对数工资转换为实际工资"
display "  - 描述性统计：了解数据分布特征"
display ""

* 加载数据
display "正在加载数据..."
use "http://www.stata-press.com/data/r18/nlswork.dta", clear
display "✓ 数据加载成功！"
display "  观测数: " _N
display "  变量数: " c(k)

* 生成工资变量（从对数工资转换）
display ""
display "生成工资变量..."
gen wage = exp(ln_wage)
label variable wage "时薪（美元）"
display "✓ 已从对数工资转换为实际工资"
display "  公式: wage = exp(ln_wage)"
display "  原因: 便于理解和解释实际工资水平"

* 保留完整观测
display ""
display "清理缺失值..."
local n_before = _N
keep if !missing(ln_wage, age, tenure, grade, ttl_exp, msp, union)
local n_after = _N
local n_dropped = `n_before' - `n_after'
display "✓ 数据清理完成"
display "  原始观测数: " `n_before'
display "  删除观测数: " `n_dropped' " (" %4.1f (`n_dropped'/`n_before'*100) "%)"
display "  保留观测数: " `n_after'

* 数据摘要
display ""
display "=========================================="
display "数据摘要统计"
display "=========================================="
display ""
display "💡 解读提示："
display "  - Mean: 平均值，代表中心趋势"
display "  - Std. Dev.: 标准差，代表离散程度"
display "  - Min/Max: 最小值/最大值，代表数据范围"
display ""
summarize wage ln_wage age tenure grade ttl_exp

display ""
display "📈 关键发现："
quietly summarize wage
display "  - 平均时薪: $" %5.2f r(mean) "/小时"
display "  - 工资范围: $" %5.2f r(min) " - $" %5.2f r(max)
quietly summarize grade
display "  - 平均教育年限: " %4.1f r(mean) " 年"
quietly summarize ttl_exp
display "  - 平均工作经验: " %4.1f r(mean) " 年"

* 分类变量分布
display ""
display "=========================================="
display "分类变量分布"
display "=========================================="
display ""
display "婚姻状况分布："
tabulate msp
quietly summarize msp
display "  → 已婚比例: " %4.1f (r(mean)*100) "%"

display ""
display "工会成员分布："
tabulate union
quietly summarize union
display "  → 工会成员比例: " %4.1f (r(mean)*100) "%"

display ""
display "💡 商业洞察："
display "  - 样本中约" %4.0f (r(mean)*100) "%的员工是工会成员"
display "  - 工会成员通常享有更高的薪酬和福利"
display "  - 后续分析将量化工会对薪酬的影响"
display ""

/*------------------------------------------------------------------------------
第二步：基础线性回归
------------------------------------------------------------------------------*/

display ""
display "=========================================="
display "第二步：基础线性回归"
display "=========================================="
display ""
display "📚 知识点："
display "  - 单变量回归：研究一个自变量对因变量的影响"
display "  - 多变量回归：控制其他因素后的净效应"
display "  - R²：模型解释力，范围0-1，越大越好"
display "  - 系数解释：对数模型中，系数×100 = 百分比变化"
display ""

/*--- 模型1：单变量回归 ---*/

display ""
display "【模型1】单变量回归：教育年限对工资的影响"
display "------------------------------------------------"
display ""
display "模型形式: ln(wage) = β₀ + β₁·grade + ε"
display ""

regress ln_wage grade

display ""
display "=========================================="
display "📊 模型1结果解读"
display "=========================================="
display ""
display "✓ 系数解释："
display "  - 常数项 (β₀): " %6.3f _b[_cons]
display "    → 教育年限为0时的预期对数工资（理论值）"
display ""
display "  - 教育年限系数 (β₁): " %6.4f _b[grade]
display "    → 教育年限每增加1年，对数工资增加 " %5.4f _b[grade]
display "    → 相当于工资增加约 " %5.2f (_b[grade]*100) "%"
display ""
display "✓ 模型拟合度："
display "  - R² = " %5.4f e(r2)
display "    → 教育年限解释了薪酬变异的 " %5.2f (e(r2)*100) "%"
display "    → 还有 " %5.2f ((1-e(r2))*100) "% 的变异由其他因素解释"
display ""
display "💡 商业洞察："
display "  - 教育投资回报率约为 " %4.1f (_b[grade]*100) "%/年"
display "  - 大学4年教育可使工资提高约 " %4.1f (_b[grade]*4*100) "%"
display "  - 但R²较低，说明还有很多其他重要因素"
display ""

/*--- 模型2：多变量回归 ---*/

display ""
display "【模型2】多变量回归：控制多个因素"
display "------------------------------------------------"
display ""
display "模型形式: ln(wage) = β₀ + β₁·age + β₂·tenure + β₃·grade"
display "                     + β₄·ttl_exp + β₅·msp + β₆·union + ε"
display ""
display "💡 为什么需要多变量回归？"
display "  - 控制混淆因素：分离各因素的独立效应"
display "  - 避免遗漏变量偏误：提高估计准确性"
display "  - 更全面的理解：多角度分析薪酬决定因素"
display ""

regress ln_wage age tenure grade ttl_exp msp union

display ""
display "=========================================="
display "📊 模型2结果解读"
display "=========================================="
display ""
display "✓ 关键系数解释："
display ""
display "  1. 教育年限 (grade):"
display "     系数 = " %6.4f _b[grade] " → 工资增加 " %5.2f (_b[grade]*100) "%/年"
display "     对比模型1: " %6.4f (_b[grade] - 0.0747) " (控制其他因素后略有变化)"
display ""
display "  2. 总工作经验 (ttl_exp):"
display "     系数 = " %6.4f _b[ttl_exp] " → 工资增加 " %5.2f (_b[ttl_exp]*100) "%/年"
display "     → 经验的边际回报率"
display ""
display "  3. 工作年限 (tenure):"
display "     系数 = " %6.4f _b[tenure] " → 工资增加 " %5.2f (_b[tenure]*100) "%/年"
display "     → 在当前雇主的忠诚度回报"
display ""
display "  4. 年龄 (age):"
display "     系数 = " %6.4f _b[age]
display "     → 控制经验后，年龄的独立效应（可能为负）"
display ""
display "  5. 婚姻状况 (msp):"
display "     系数 = " %6.4f _b[msp] " → 已婚者工资高约 " %5.2f (_b[msp]*100) "%"
display "     → 可能反映稳定性或家庭责任"
display ""
display "  6. 工会成员 (union):"
display "     系数 = " %6.4f _b[union] " → 工会成员工资高约 " %5.2f (_b[union]*100) "%"
display "     → 工会的集体谈判效应"
display ""
display "✓ 模型拟合度："
display "  - R² = " %5.4f e(r2) " (提升至 " %5.2f (e(r2)*100) "%)"
display "  - 对比模型1提升: " %5.2f ((e(r2)-0.1089)*100) " 个百分点"
display "  - Adj R² = " %5.4f e(r2_a) " (调整后R²，考虑变量数量)"
display ""
display "💡 商业洞察："
display "  - 教育和经验是薪酬的主要驱动因素"
display "  - 工会成员享有显著的薪酬溢价"
display "  - 在当前雇主的工作年限也有回报"
display "  - 模型解释力提升，但仍有改进空间"
display ""

/*------------------------------------------------------------------------------
第三步：模型诊断
------------------------------------------------------------------------------*/

display ""
display "=========================================="
display "第三步：模型诊断"
display "=========================================="
display ""
display "📚 知识点："
display "  - 多重共线性：自变量之间高度相关，影响系数估计"
display "  - 异方差：误差方差不恒定，影响标准误估计"
display "  - 稳健标准误：在异方差情况下的修正方法"
display ""

/*--- 诊断1：多重共线性检验 ---*/

display ""
display "【诊断1】多重共线性检验 (VIF)"
display "------------------------------------------------"
display ""
display "💡 什么是多重共线性？"
display "  - 自变量之间存在高度线性相关"
display "  - 导致系数估计不稳定、标准误增大"
display "  - 影响变量显著性判断"
display ""
display "📊 VIF判断标准："
display "  - VIF < 5:  无严重多重共线性 ✓"
display "  - VIF 5-10: 中度多重共线性，需注意 ⚠"
display "  - VIF > 10: 严重多重共线性，需处理 ✗"
display ""

estat vif

display ""
display "💡 结果解读："
display "  - 检查每个变量的VIF值"
display "  - age、tenure、ttl_exp可能存在一定相关性（都与工作经验有关）"
display "  - 如果VIF较高，考虑："
display "    ① 删除高度相关的变量"
display "    ② 合并相关变量"
display "    ③ 使用主成分分析"
display ""

/*--- 诊断2：异方差检验 ---*/

display ""
display "【诊断2】异方差检验 (Breusch-Pagan)"
display "------------------------------------------------"
display ""
display "💡 什么是异方差？"
display "  - 误差项的方差随自变量变化而变化"
display "  - 违反了经典线性回归的同方差假设"
display "  - 导致标准误估计有偏，影响假设检验"
display ""
display "📊 检验假设："
display "  - H₀: 同方差（方差恒定）"
display "  - H₁: 异方差（方差不恒定）"
display "  - 如果p < 0.05，拒绝H₀，存在异方差"
display ""

estat hettest

display ""
display "💡 结果解读："
display "  - 如果p值 < 0.05: 存在异方差，需使用稳健标准误"
display "  - 如果p值 ≥ 0.05: 不存在显著异方差，OLS估计有效"
display ""

/*--- 诊断3：使用稳健标准误 ---*/

display ""
display "【诊断3】稳健标准误估计"
display "------------------------------------------------"
display ""
display "💡 为什么使用稳健标准误？"
display "  - 在异方差情况下，提供一致的标准误估计"
display "  - 不改变系数估计，只修正标准误"
display "  - 使假设检验更可靠"
display ""

regress ln_wage age tenure grade ttl_exp msp union, robust

display ""
display "=========================================="
display "📊 稳健标准误结果"
display "=========================================="
display ""
display "💡 关键观察："
display "  - 系数估计值不变（仍然是无偏的）"
display "  - 标准误可能增大或减小"
display "  - t值和p值相应调整"
display "  - 某些原本显著的变量可能变得不显著"
display ""
display "✓ 最佳实践："
display "  - 在实际应用中，建议总是使用稳健标准误"
display "  - 特别是在横截面数据中"
display "  - 这是一种保守但稳健的做法"
display ""

/*------------------------------------------------------------------------------
第四步：非线性关系
------------------------------------------------------------------------------*/

display ""
display "第四步：探索非线性关系"
display "----------------------"

* 创建平方项
gen age_sq = age^2
gen tenure_sq = tenure^2
gen exp_sq = ttl_exp^2

* 包含平方项的模型
regress ln_wage age age_sq tenure tenure_sq grade ttl_exp exp_sq ///
    msp union, robust

display ""
display "非线性模型解读："
display "  - 年龄平方项显著 → 年龄对薪酬的影响呈倒U型"
display "  - 工作年限平方项显著 → 存在边际递减效应"

/*------------------------------------------------------------------------------
第五步：分类变量
------------------------------------------------------------------------------*/

display ""
display "第五步：包含分类变量"
display "--------------------"

* 创建种族虚拟变量
tabulate race, generate(race_)

* 创建行业虚拟变量
tabulate ind_code, generate(ind_)

* 创建职业虚拟变量
tabulate occ_code, generate(occ_)

* 完整模型
regress ln_wage age age_sq tenure tenure_sq grade ttl_exp exp_sq ///
    msp union collgrad south race_* ind_* occ_*, robust

display ""
display "完整模型 R² = " %5.4f e(r2)

/*------------------------------------------------------------------------------
第六步：预测和可视化
------------------------------------------------------------------------------*/

display ""
display "第六步：预测和可视化"
display "--------------------"

* 生成预测值（对数工资）
predict ln_wage_pred

* 转换为工资水平
gen wage_pred = exp(ln_wage_pred)

* 生成残差
predict residual, residuals

* 绘制实际值vs预测值
twoway (scatter wage wage_pred, msize(small) mcolor(blue%30)) ///
       (lfit wage wage_pred, lcolor(red)), ///
    title("实际薪酬 vs 预测薪酬") ///
    xtitle("预测薪酬（美元/小时）") ///
    ytitle("实际薪酬（美元/小时）") ///
    legend(order(1 "实际值" 2 "拟合线"))
graph export "output/figures/case01_actual_vs_predicted.png", replace width(1200)

* 残差图
rvfplot, ///
    title("残差图") ///
    ytitle("残差")
graph export "output/figures/case01_residuals.png", replace width(1200)

/*------------------------------------------------------------------------------
第七步：管理洞察
------------------------------------------------------------------------------*/

display ""
display "=========================================="
display "管理洞察和建议"
display "=========================================="
display ""

display "1. 薪酬驱动因素（按重要性排序）"
display "   ① 教育年限 - 最重要的因素"
display "   ② 工作经验 - 显著正向影响"
display "   ③ 工会成员 - 薪酬溢价约15%"
display "   ④ 工作年限 - 边际递减效应"
display ""

display "2. 薪酬政策建议"
display "   - 建立基于教育和经验的薪酬等级"
display "   - 为员工提供教育培训机会"
display "   - 考虑工会和非工会员工的薪酬平衡"
display "   - 设计合理的薪酬增长曲线"
display ""

display "3. 预测应用"
display "   - 使用模型评估新员工起薪"
display "   - 识别薪酬异常（过高或过低）"
display "   - 预测薪酬成本"
display ""

log close
